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abstract 

We study the thermally activated oscillations, or capillary waves, of a neutral metal cluster within 
the liquid drop model. These deformations correspond to a surface roughness which we characterize 
by a single parameter A. We derive a simple analytic approximate expression determining A as a 
function of temperature and cluster size. We then estimate the induced effects on shell structure by 
means of a periodic orbit analysis and compare with recent data for shell energy of sodium clusters in 
the size range 50 < < 250. A small surface roughness A ~ 0.6 A is seen to give a reasonable account 
of the decrease of amplitude of the shell structure observed in experiment. Moreover - contrary to 
usual Jahn- Teller type of deformations - roughness correctly reproduces the shape of the shell energy 
in the domain of sizes considered in experiment. 
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1 Introduction 



Since the discovery of shell effects in metal clusters, the mean field approach with delocalized electrons 
has been a very efficient tool for describing a wide variety of phenomena : shell and supershell effects, 
dipole polarizability and optical excitations, fission . . . (for reviews see Refs. 

This type of approach mainly relies on the jellium model where the ionic background is considered 
as a smooth and uniform distribution of positive charges. It is best legitimate for simple metals 
(and to a lesser extend for noble metals) with delocalized valence electrons, almost insensitive to the 
actual arrangement of the ionic cores. Hence the best candidates for this approximation are the alkali 
metals, as can be inferred from the quasi-sphericity of their Fermi surface (revealing a weak interaction 
between ionic cores and valence electrons) . As a result the compressibility of these solids is close to its 
electron gaz value, and the surface tension is correctly described by the jellium models (the agreement 
being better for small electronic density, see Ref. Q). 

It was realized early that deformation effects had to be taken into account for a realistic des- 
cription of metal clusters (see eg. the review Q]). Within density functional theories, this can be 
achieved by imposing to the jellium the shape that suits the electrons best [5-^|. Deformations have 
also been studied in less elaborate models (deformed external mean-field) [|-13] and there is a good 



overall agreement with experimental data for ionization potentials, dissociation energies and splitting 
of dipole resonances for relatively small clusters (less than 40 atoms) . 

The above mentioned deformations are of Jahn- Teller type and occur between major shell closures, 
where lowering the symmetry leads to a gain in energy. Another type of surface deformation has also to 
be considered which consists in surface irregularities of very large multipolarities. These deformations 
not only lower the shell effect but also introduce randomness into the spectrum. This was first noticed 
by Gor'kov and Eliashberg [I4| who claim that "the distribution of the levels should be random even 



if the particles have the same volume and a good shape, say spherical particles of equal size. The 
point is that electrons in the metal have a wavelength of the order of atomic dimensions. Therefore 
surface irregularities of atomic size are sufficient to make the level distribution perfectly random". 
This statement has to be tempered in view of the success of the jellium model. Nevertheless it makes 
no doubt that the surface of a cluster has atomic size irregularities, it is important to estimate their 
amplitude and to evaluate the resulting effect on the physical observables. 

The fact that disorder is located on the surface is legitimated because the elastic mean free path 
of an electron in the bulk is typically of order of several hundreds of Angstroms, whereas an electron 
experience collisions on the surface of the cluster about each 10 A. Indeed it was shown in Ref. [^] 
that the scattering of electrons on the fluctuation of the positive ions had an effect of the same order 
as that of thermal distribution of occupancy probability ; which in turn is shown in the present work 
to be negligeable compared to the effect of shape fluctuations. More microscopically, bulk disorder 
would be represented by fluctuations of the bottom of the potential well and in the large size limit high 
lying states tend to be insensitive to this perturbation, whereas the effects of surface disorder increase 



when one goes up in the spectrum |16|. Note that such irregularities are to be taken into account 
also when the cluster is "liquid-like" : the mean velocity of the ionic cores is always by several orders 
of magnitude smaller than the typical electronic Fermi velocity. Hence, as far as electronic motion is 
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concerned, the ionic cores can be considered as frozen and this necessarily imphes a certain degree of 
surface roughness. 

In the present paper we use a hquid drop model to study the thermally activated surface defor- 
mations, or capillary waves, of a neutral spherical cluster (Sec. 2). These deformations correspond 
to a surface roughness which we characterize by a single parameter A. We derive a simple analytic 
approximate expression determining the behaviour of A as a function of temperature and cluster size. 
Then in Section 3 we discuss the influence of shape fluctuations on the level density using a trace 
formula in rough billiards and compare with thermal effects linked with Fermi occupation number 
of the energy levels (for this purpose we give the general form of shell corrections in the presence of 
temperature in Appendix C). Finally we present our conclusions and discuss possible refinements of 
our approach in Sec. 4. 



2 Liquid drop model 



In the liquid drop model, a cluster is described as a droplet of incompressible fluid whose shape can 
be parameterized by a set of normal coordinates a^/x obtained by expanding the surface in spherical 
harmonics [171,181: 



r{n,t) = R 



1 + 



(1) 



The r.h.s. of Eq. (0) is made real by imposing a\-^ = i—)^o:\^. The summation stops at a 
Debye cutoff A estimated by equating the number of surface modes to the number of atoms on the 
surface, this yields A = (3\/47rA^)^/^ ~ 2.20 N^^^. The droplet being considered as incompressible 
one should impose volume conservation. If the cluster contains N atoms it should have a volume 
V = 47rR^/3 with R = rs N^^^ (rg being the Wigner-Seitz radius of the material). This leads to the 
relation aooV47r = — J2xfi |q^A/xPi valid to leading order. Also the modes A = 1 which correspond to a 
global translation of the drop should be omitted in the summation (|^). 

Eq. (|l|) yields a kinetic energy T and a surface energy I4urf = crA, where a is the surface tension 
and A the surface area corresponding to (||). Including terms up to second order in the a's one 



obtains |13,|T|1 : 
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and 



surf = 47ri?V + ^ E I"AmI'(A - 1)(A + 2) 



(2) 



where po is the specific mass of the material considered. 

One can also take into account a curvature term in the potential energy 

In (^) 7 is an intrinsic curvature energy parameter, TZi and 7^2 are the principal radii of curvature. 
It turns out (see Appendix A) that taking this term into account exactly amounts to replacing in Eq. 
(^) the surface tension by an effective term a — > a* = a + j / {2R) with which we will work henceforth. 



(3) 
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No other contribution to the potential energy has to be taken into account because we consider 
concomitant deformations of the jehium and of the valence electron cloud of a neutral cluster (hence 
there is no other electrostatic deformation energy than the one included in (|2|,^)). We neglect at 
this level finite size quantum effects. We use for the surface tension a the value of the bulk material 
extracted from experiment in Ref. and this implicitly contains quantal effects associated with the 
kinetic energy of the electrons near the surface. Hence including quantum effects in the present 
description would double count this contribution ; the appropriate procedure would be to use a 



Strutinsky shell correction, cf. the discussion at the end of the paper. In Ref. |20| a description 
analogous to the present one has been shown to account accurately of the monovacancy formation 
energy in simple metals such as the one we are interested in. This gives us confidence in the ability of a 



liquid drop model to describe atomic size irregularities. Note that in |2C] the value of a is renormalized 
in order to describe an ideally flat surface. This procedure should not be employed here because we 
want the surface tension of a large cluster to tend to the one of the bulk material. 

Eqs. dH) and (^) correspond to a liquid drop Lagrangian £^13 [a ,ct\^] = T — Kiurf — V'curv with 
normal modes ax^{t) = Ax^exp{iu;xt) of pulsation uj\ given by : 

u;2 = A(A-l)(A + 2)^, (4) 

and the classical energy of the mode is Ex^ = a*B?{\ — 1)(A + 2)|74A/ip. The average value of 
the amplitude |^A/^P of the thermally activated mode is determined by writing Exf_i = ksT. We 
use here classical statistical mechanics, the quantal analog would be Ex^ = huJx{nx + 1/2), where 
nx = [ex.p{hu}x/kBT) — l\~^ is a Bose occupation factor. Such a description has been used for describing 



the surface oscillations of liquid helium |21|, but here the motion of the surface is classical : kgT ^ huj\ 
(from (^) hLJ\ ~ 130 K for sodium). 

From (HI) the quantity r(J7,t) averaged over the surface has a mean value R{1 + aoo/\/47r) and a 
standard deviation A which is is given by = R^J2x>2 l'^A^iP/(4'/r), hence A is independent of time. 
The explicit formula reads : 

^ 2A + 1 ksT (2A-l)(2A + 5) 

A — z 

In the r.h.s. of Eq. (|5|) we replaced the discrete summation by the first term of its Euler-MacLaurin 
expansion. Figure (la) displays the result of Eq. (^) for sodium clusters at temperatures T = 200 K 
and T = 450 K in the size region 20 < < 1000. It is difficult to determine the precise value of a* 
to be used in Eq. (|5|) : the surface and curvature parameters a and 7 depend on temperature and 
of the actual phase (liquid or solid) of the aggregate. Hence we used several values of a and 7 : a 
lower bound for A is obtained by taking the values a = 190 K.A~^ (which is the solid-vapor value 
extrapolated to zero temperature in Ref. |9|) and 7 = 285 K.A'^ [|o|. The upper bound is obtained 
by taking 7 = and a = 145 K.A~^ (which is the liquid- vapor surface tension at melting |]l^). These 
values of a correspond to a droplet parameter = 47rr|(T which ranges from 0.68 eV (for a = 145 
K.A-2) to 0.89 eV (for a = 190 K.A^^). Indeed one can find a large dispersion of Ug in the literature : 
the value 0.54 eV was used in Refs. [l^,^ (from a fit to theoretical values of clusters' energy) ; in 
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Ref. 1^ the value = 0.7 eV was extracted from the bulk surface tension and in Ref. |23] the value 
as = 1.02 eV was obtained via experimental determination of clusters' cohesive energy. 



(a) (b) 




N N 



Figure 1: Part (a) : A as a function of N for sodium clusters at temperatures 200 K and 450 K, 
from Eq. (I). The shaded zones correspond to different values of the surface tension and curvature 
parameter (see the text). Part (b) : same as part (a) withdrawing in Eq. (||) the contribution of the 
two first multipolarities (A = 2 and A = 3). 

As several studies in the field have demonstrated (see e.g the work of Yannouleas and Landman 
[p^ , ) the liquid-drop Lagrangian C^o is not adequate for describing the deformations of small 
multipolarities, which are determined mainly by shell effects (see the discussion at the end of the 
paper). In order to estimate the role of these multipolarities in the amplitude of the roughness we 
have redrawn Figure (la) withdrawing the contributions of A = 2 and 3 in the summation (^). The 
result is displayed on part (b) of Fig. 1. As one should expect this significantly reduces the value of 
A. For instance, at T = 450 K the typical roughness for ~ 200 is reduced from 1 A to 0.85 A. 

From the discussion above and the comparison of Figures (la) and (lb), we see that due to the 
simplicity of our model we cannot precisely determine the surface roughness of sodium clusters at 
typical experimental temperatures. For sodium rg = 2.08 A (thermal variations of rg are negligeable 
here), and we can only say that typical roughnesses are of order of 30 to 50 % of the interatomic 
distance. However it makes no doubt that there is a thermal activation of capillary waves which gives 
a contribution to the surface roughness of the type given by Eq. (^) for large enough multipolarities. 
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We will see below that this has an important effect on shell structure, and that this effect is crucial 
for understanding the shape of the shell energy determined in experiment. 

Note that an estimation of the geometrical corrugation of a solid surface at zero temperature [^] 
yields values of A of order of 20 % of rg. Hence the reduction of shell oscillations presented below is 
a general phenomenon which does not depend on the solid or liquid structure of the aggregate. 



3 shell energy in a rough sphere 



The influence of surface roughness on the level statistics has been discussed in Refs. pq-pq] and more 
recently in a 2D model |^9| . In the present work we concentrate on its effect on level density and shell 
structure. Similar effects have been recently studied in Refs. [30, 31| and The spirit of the present 
section is very similar to the one of Ref. which presents numerical results in a corrugated mean 
field. However the qualitative conclusions are different : in Ref. corrugation is seen to imply a shift 
in the supershell structure. This effect - not seen in the present study - seems to be due to the fact 



that, in a finite depth potential such as the one used in [32|, roughness leads to an effective mean-field 
where the phase difference between the orbits is modified. In any case, for the small roughnesses we 



are using here, the shift in the shell structure found in Ref. |32] is small. In Ref. |30| the disorder is 
modelized by the addition of a random Hamiltonian to the mean field and the results are comparable 
to the one presented below. The approach has been further extended in Ref. where the effects of 
disorder on energetics of lithium, sodium and potassium clusters where taken into account in a liquid 
drop plus shell-correction model. Here the Hamiltonian for the deformation is only of liquid drop 
type ; however we do take into account the thermally activated oscillations for the Hamiltonian we 
consider (see the discussion in Sec. 4). 



In the present work the electrons are considered to be moving independently in an infinite 
potential well (a billiard) having a shape approximatively spherical (as given by dl])). Hence we can 
consider that the actual shape is obtained by a random deviation from a perfect sphere, with Gaussian 
fluctuations of standard deviation A determined above. The choice of Gaussian fluctuations reflects 
the fact that the distribution of each is Gaussian (according to classical mechanics) . Then invoking 
the central limit theorem it can reasonably be considered that the shape fluctuations are of Gaussian 
type. We consider an ensemble of clusters (such as one would expect in a molecular beam) and we 
will present results averaged over this ensemble. The radius R of the average sphere scales with so 
that the mean electronic density is kept constant and equal to its bulk value : R = Vg N^^^ . 

The level density in a rough billiard with small size surface irregularities was studied in Ref. [^] 
and a semiclassical trace formula averaged over surface disorder was derived. The important feature 
of the level density is the gradual disappearance of shell effects with increasing energy : near the Fermi 
level the electronic wavelength is of order of the typical size of the surface defects and the induced shift 
of the eigenlevels leads after averaging to a structureless level density. The bottom of the spectrum 
is not affected because low lying state have a wavelength much larger than the surface perturbations 
(accordingly, the effect on level statistics is different at the bottom of the spectrum and near the Fermi 
energy [||]). 

It is shown in Ref. that the oscillatory part of the electronic energy i?sheii (the so-called shell 
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energy) can be expressed on average as a sum over classical periodic orbits in a perfect sphere 



^shcii(iV, A) ~ 



2m 



sin(A:^L + utt/2) exp{-2n(A:^A)^ cos'-" 6} 



PO 



(6) 



In (^) m is the electron mass, kp is the smooth part {i.e. non-oscillatory) of the Fermi wave-vector, 
which is to a good approximation equal to the bulk wave- vector Kp = r^^(97r/4)^/'^. Egheii in ® is a 
quantity averaged over surface disorder, but the sum is performed over all the periodic orbits (POs) 
of a perfect sphere (see [16|). L is the length of a PO, A is an amplitude slowly depending on kp, v 
is a Maslov index, n is the number of bounces of the PO on the sphere and 9 is the bouncing angle ; 
all these quantities depend on the PO considered, see Appendix B for further details. When A = 0, 
@ follows from Balian and Bloch's trace formula for the sphere Since kp is nearly constant, the 
main A^-dependence in (^) is due to the scaling of the cluster's size according to = rg A^^/^ (L scales 
like R,A(x R^/"^ or for some orbits, see Appendix A). 
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Figure 2: £'sheii(A, A) as a function of N in sodium clusters. The lower plot compares the experimental 
results of Ref. |^4j (black points) with the values obtained by fixing A = 0.9 A (dashed line) and 
A = 0.63 A (solid line). The upper plot compares the experimental data with the values obtained by 
determining for each cluster size A via (||), taking a = 190 K.A~^ and 7 = 285 K.A~^, withdrawing 
the contributions of A = 2 and 3. The dashed line corresponds to a temperature T = 400 K and the 
sohd hue to T = 300 K. 
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-Bsheii as given by (0) is compared on Fig. 2 with the value determined by Chandezon et al. in 



Ref. [34|. In this reference an evaporation model is used for extracting the shell energy from the 
abundance distribution in a beam and results for clusters of sizes ranging from 50 to 230 atoms are 
obtained. This experiment is important for our study because it concerns relatively large clusters 
and our approach is limited to this domain for the two following reasons : (i) we use a semiclassical 
approach more accurate for large sizes (see for instance the comparison with exact results on Fig. 4) 
and moreover (ii) the macroscopic concept of roughness is meaningless for small sizes ; for instance 
in a cluster with = 20 atoms the Debye cut-off fixes the maximum angular momentum of surface 
deformation to be A ~ 2.2 N^/^ ~ 6 and in this regime the concept of roughness is of marginal 
importance (see below the discussion of Jahn- Teller effect). In the following we state rather loosely 
that our approach is relevant for sizes > 100. 

On the lower part of Fig. 2 the dashed curve corresponds to a constant value A = 0.9 A (inde- 
pendent of A'^) which - from Fig. 1 - is a typical value at T = 450 K (this temperature is in agreement 



with usual evaporation conditions [35-37|). The fact that this value of A leads to a too large damping 
of shell structure should not worry us at this level : as stated in the previous section, the liquid drop 
model does not accurately determine the value of A because it does not properly describe deformations 
of small multipolarities. These deformations should be described with a more elaborate procedure and 
may be less easily thermally excited (see the discussion of Sec. 4). This is confirmed by the solid curved 
which is drawn for A = 0.63 A and which gives a better account of the data. Note the sensitivity of 
-E'sheii to the value of A/rg : formula (P) has a schematic large N behaviour of the form : 

i?sheii(iV, A) ~ e^ivVe exp(-AVr|) sin{N^I^) . (7) 

where £p = h'^Kp/{2m) is the bulk Fermi energy. For clarity we have dropped in the sine and the 
exponent of ^ important but dimensionless factors. This will be done also in ; the derivation of 
these formulae is explained in Appendix C (Eq. (|C9|)). Hence the shell structure is very sensitive to a 
small surface roughness ; this point and the validity of formulae of type @ have been further tested 



on a numerical example in Ref. |38]. 



The value A = 0.63 A does not quite correspond to the estimation of Fig. 1 for T = 450 K : 
as explained in the conclusion the liquid drop Lagrangian seems to underestimate the stiffness of the 
potential for the deformation parameters. It may also happen that the electrons experience a mean- 
field which - due to the diffuseness of its surface - is less corrugated than the ionic background. In the 
same line, instead of fixing A to a constant value, one should according to (|5|) take size-dependence 
of the roughness into account. This improves the agreement with experiment for ~ 50 since in 
this region the value of A decreases significantly (see Fig. 1) and this leads to lower damping of the 
theoretical curve which comes closer to experiment. This has been done on the upper part of Fig. 2 
which is drawn in the case a = 190 K.A~^ and 7 = 285 K.A~^ ; these values have been chosen because 
they lead to small values of A and to relatively good agreement with experiment. The dashed line 
correspond to T = 400 K and the solid line to T = 300 K. For each temperature and cluster size A 
was determined via (^), withdrawing the contributions of A = 2 and 3. However, such a refinement 
is unnecessary for larger cluster sizes in view of the small A^-dependence of (^) for large A^. Besides, 
the simple model with a constant value of A = 0.63 A gives already a satisfactory agreement with 
experimental data. 
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In fact our approach is more strongly supported by the very good comparison of theory an ex- 
periment for the shape of the curve for the shell energy than for the comparison with the amplitude. 
Indeed the agreement with the amplitude may not be as good as presented on Fig. 2 because there 
should be some room left for an extra reduction of the amplitude of the shell energy due do quan- 
tum mechanically driven deformations, corresponding to relatively small multipolarities (A = 2 or 3 
typically). Nevertheless, due to the high sensitivity of formula (^) to small changes of A, we still can 
conclude from Fig. 2 that the typical roughness is of order of 0.6 A. 



Concerning the shape of the curve, the experimental results of Ref. |34] are surprising because 
they are in contradiction with the common belief that cluster's deformations are only governed by 
Jahn- Teller effects. For instance, the dissociation energies and ionization potentials of simple metal 
clusters of relatively small sizes are well accounted for by models where the Jahn- Teller effect is the 



only mechanism of deformation [12,2^ (the agreement with experiment survives up to size N ~ 100 
for the ionization potentials of potassium, see Ref. |^^). This phenomenon was expected to occur even 
for large values of N, see e.g. the zero temperature results of Refs. [p|-pll,[l3|| , or the finite temperature 



results of Ref. [40|. In these studies, the deformations occur between shell closure and their main 
effect is to remove the upper part of the shell oscillations ; the shell energy is predicted to have sharp 
negative spikes in vicinity of the magic numbers (these spikes correspond semiclassically to long POs). 
On the other hand surface roughness suppresses long POs and reduces shell structure more uniformly, 
as seen in the experimental data of Ref. Hence we feel that previous theoretical approaches 

overestimate the role of the Jahn- Teller mechanism : the very specific shape of shell energy they 
predict is not seen in the experiment of Chandezon et al. The comparison between our approach and 
the experimental results for the sizes > 100 firmly establishes that there is a qualitatively important 
effect of roughness. 

On the quantitative level, one can also notice that typical theoretical studies overestimate the shell 



effect for clusters of large sizes : compare Fig. 3 of Ref. |3J] with similar figures of Refs. |^11, 
Temperature effects improve the agreement (see Ref. [^]) but there is still a mismatch of order of 40 
% for the amplitude of shell oscillations (see Fig. 3 of Ref. [Q), leaving room for improvement due 
to surface roughness. 

For further comparison with experimental data we display on Table 1 the magic numbers in the 
region N < 1300. The first column shows the A = results from the semiclassical formula The 
second column displays the exact result in the perfect sphere and merely tests the accuracy of the 
semi-classical periodic orbit expansion used in the first column. Note that the magic numbers of these 



two columns are almost identical to the results of Bulgac and Lewenkopf |1C| who use a quadrupole 
deformation of a spherical billiard model within the shell correction method : this is due to the fact 
that, as stated above, there is no Jahn- Teller deformation at shell closure. In the third column we 
show the magic numbers obtained from (^) with A = 0.63 A. The three first columns compare well 
with the experimental ones from Chandezon et al. (column 4), and roughness has only a small effect 
on the location of the magic numbers. We still produce these data because they justify the billiard 
model we are using : the magic numbers from Table 1 are in better agreement with experiment than 
the one obtained with harmonic oscillators |^,11| or more elaborate potentials |13, 40|. Hence, as far 
as the phase difference between the contribution of its POs to (^) is concerned, the billiard model is 
presumably close to the experimental situation since it allows a good prediction of the minima in the 
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Table 1: Magic numbers in the perfectly spherical billiard (column 1 : PO expansion, column 2 : exact 
results) and in the rough billiard (column 3). Column 4 shows the experimental results of Chandezon 
et al. 

shell energy. However the electrons are sensitive to a mean-field which - due to the finite range of the 
electron/ion and electron/electron interaction - could be less corrugated than the ionic background. 
This would lead to an effective decrease of A and may help improving the model by reducing the 
importance of the ionic corrugation, leaving some room for an extra decrease in amplitude due to 
deformations of Jahn- Teller type. 

Note that in the present treatment the effects of temperature are indirect : although the usual 
temperatures reached in experiments are small compared to the Fermi energy (one remains in the 
highly degenerate limit fegT <C Sp) they are sufficient to induce a disorder of the cluster's shape 
which has a sizeable effect on shell structure. For comparison one can derive a formula (similar to (^)) 
encompassing the effect of a Fermi occupation function in the energy levels of the electron gas. The free 
energy F{N, T) is more appropriate than the total energy for evaluating these effects. Indeed, based on 
Weisskopf 's approach, the electronic contribution to the evaporation rate of a neutral monomer from 
a cluster of size N can be shown to be approximatively proportional to exp{[F(A^) — F{N — 1)]//cbT} 
[^,^. The general formula for the oscillating part of the free energy is derived in Appendix C and 
reads in the case of a billiard : 



i^shell(iV,T) ^^^^i^sin(fc^L + i.vr/2)Fi(X) , (8) 
PO A* 
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where A;^ is the non-oscillatory part of the quantity fc^ defined by = h^kf^/lm, fi being the che- 
mical potential. Again is to a good approximation equal to the bulk Fermi wave vector Kp. 
X = {-K /2)tLkp /kfj, is a dimensionless quantity which tends to zero at T = (r = ksTjep is the 
reduce temperature). More precisely it can be considered as small if the thermal wave length At = 
{271^? / mkgTY /'^ is large compared to y/rgL {X = 2'K'^L/k^\j). Fi is a dimensionless damping function 
defined in Appendix C (Eq. Q). 

We compare on Fig. 3 the effects of a temperature T = 750 K, with those of a constant roughness 
A = 0.63 A. The shell energy is displayed as a function of N^^^ for sodium clusters of size < 
3400. There is a striking difference with the A/"-dependence obtained via usual temperature effects on 
occupation numbers. As one notices from the figure, roughness damps the oscillations in the total 
energy with an overall factor of the type exp(— A^/r|) (c/ (0)) without modifying the qualitative 
features of the supershells ; whereas temperature leads to a A^-dependent damping of schematic form 
(c/ Appendix C, Eq. @) : 

Fsheii(Ar,r) ~ epN'/^Fi{TN'/^)sm{N^/') , (9) 

(as in Eq. (0) we have omitted numerical factors in the sine and Fi function). Hence the effects of 
thermal distribution of occupation numbers is to wash out the beating pattern of the shell energy by 
exponentially damping the large N oscillations (see also Fig. 4). 
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Figure 3: E'sheii (expressed in eV) as a function of N^^^ in sodium clusters. The upper plot is obtained 
by taking a constant roughness A = 0.63 A. The lower plot corresponds to formula ( |C^ ) for a perfectly 
spherical aggregate with a temperature r = kgT/ep = 0.02, i.e. T = 750 K. 
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Note also the efficiency of a small roughness for diminishing shell effect : for clusters of size 
800 < N < 1200, from Eqs. (^) and (^), one can see for instance that the effect of a small roughness 
A = 0.2 rs ^ R/50 on i^sheii is similar to the one of a temperature of about 550 K on Fsheii- In 
the range 300 < N < 500 the same roughness corresponds to a temperature T ~ 750 K, and for 
50 < iV < 250 it corresponds to T ~ 1000 K (around = 200) or 1400 K (around N = 100). As a 
result, in a range of sizes and temperatures commonly reached in experiment (N ~ 500 and T ~ 400 
K) the effect of roughness on shell structure is dominant compared to that of thermal distribution 
of occupation numbers. Furthermore, as discussed above, it seems from the experimental results of 
Ref 1 34 1 that in about the same region, Jahn- Teller deformations of small multipolarity play a smaller 
role in the shape of the shell energy than predicted by usual theoretical studies (see Fig. 3 of Ref. p^). 
We present here roughness as a concomitant phenomenon which (according to comparison with the 
data of Ref. [34]) seems more relevant for large clusters. As discussed in next section both phenomena 
should be taken into account for a proper description of deformations of large clusters. 



4 Discussion 

One of the original interests of metal clusters was to provide a physical realization of a discrete and 
random spectrum. It was long thought that the randomness of the levels would lead to a structureless 
level density and theoretical works were mostly devoted to the study of the two point form factor of the 
spectrum [^,^. It is only during the last two decades that molecular beams have made it possible 
to work in a regime where the size of the cluster is well defined and small compared to the electronic 
mean free path. In this regime one observes shell effects as a result of finite size quantum effects [l|,^. 
Nevertheless the remark of Gor'kov and Eliashberg quoted in the introduction remains valid to some 
extend and the present work aims to reconcile these two views by providing a semi-classical description 
of shell structure in presence of shape disorder. 

The model we have considered is schematic ; a more elaborate procedure would be to design a 
Lagrangian for the surface deformation encompassing the effects of shell structure (in the spirit of 
Strutinsky shell corrections) : 

£ = CLD[a\ti , ax^] - Eshcii[a\fj] ■ (10) 

This approach is typical in the study of deformations of finite fermionic systems. It is used for 
determining the equilibrium shape, i.e. the set of minimizing the total potential energy. Between 
shell closure it leads to a ground-state in which the equilibrium value of some of the a's is non-zero 
(mainly for small A), contrary to what is obtained in Sec. 2 where all the a's are zero at equilibrium. 



Using the Lagrangian (10) would also modify the stiffness of the potential near the minimum (since a 
term -Bgheii would be added to the potential used in Sec. 2). 

Following the procedure of the present paper one should go one step further and study the thermally 
activated vibrations of the normal modes in the Hamiltonian (p^). Hence the usual Jahn- Teller 
deformations correspond to the first step of the procedure just exposed and to small multipolarity 
whereas surface roughness to the second step (and to large multipolarity). 

It would be of great interest to verify if the agreement obtained in Sec. 3 with experimental values 
would persist when describing surface oscillations with a Lagrangian such as C defined in Eq. (^). This 
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form of Lagrangian is legitimated by confrontation with experiments in the small and intermediate size 
domain, where it is commonly admitted that clusters of size < 100 experience static deformations 
of small multipolarities [^,^]. The success of the present model (which uses C^o and do not include 
such Jahn- Teller deformations) in the size range > 100 might be explained by the decrease of the 
shell-energy contribution to C due to an intrinsic roughness of the surface. Schematically one might say 
that there is less difference between a rough sphere and a rough ellipsoid than between a perfect sphere 
and a perfect ellipsoid. A similar phenomenon explains the disappearance of Jahn- Teller deformations 
with increasing temperature, see Ref. [^] where deformation is seen to be suppressed by thermal 
fluctuations. Note however that the phenomenon predicted in this reference is size-dependant, i.e. 
not uniform for all cluster sizes (as roughness would be), see the precise discussion in pO| . 



Note finally that thermally induced shape fluctuations have already been investigated in the study 
of the broadening of plasma resonances of metal clusters in Refs. |44-47| and very recently in the 
study of shell structure for clusters of size smaller than A^ = 100 The ideas are similar to the 

one exposed above, however the allowed deformations are limited to simple shapes, whereas it has 
been considered necessary in the present work to include deformations of very large multipolarities for 
studying surface roughness. 
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Appendix A 

In this appendix we compute the curvature energy Vcmv (deflned in Eq. (^) for a droplet of shape given 
by (|l|). This amounts to evaluate the integral C = J dA{l/TZi + I/IZ2) for a boundary approximatively 
spherical. This can be done by noticing that if A is the surface area of a given boundary, the 
modification 5 A caused by an infinitesimal displacement of the boundary reads 5 A = / dA5C,{l/lZi + 
1/7^2); where 6C, is the normal segment between the undeformed boundary and the deformed one (see 
eg. chap. VII). If this modification corresponds to a modification of Eq. (||) by r — > r + 5r, one 
can compute 5A and 5C, in terms of 5r. This allows to write C in the form 

where K{Q) = + (dgr)^ + (S^r)^/ sin^ 9. Then, writing r{^) = R[l + h{Q)] and neglecting terms of 
order greater than 0{h'^) one obtains : 

C = 2R [ dn !l + h+ lidehf + (d^hf] + 0{h^) . (A2) 

J L 2 2 sm y J 

The surface area can be expressed in a similar manner (see ||4^) : 
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A= fdA= f dnrK^/^ = R^ fdn ((1 + hf + lidehf + — ^ {d^A + 0{h^) . (A3) 
J J J [ 2 2 sin 6 J 

The condition of volume conservation imposes / dO, (1 + /i)^ = JdQ (1 + /i) + 0{h^). Hence, com- 
paring Eqs. ( [A2| ) and ( |A3D one sees that for small deformations the curvature integral is proportional 
to the surface area : C = 2A/R + 0{h^). The corrections are of third order in the deformation, they 
are given in Ref. (chap. 6) for spheroidal and harmonic deformations. Here we are interested only 
in terms up to order 0{h?)^ thus the curvature energy Vcurv = 7^/4 is equal up to a multiplicative 
constant to the surface term I^urf = a A : V^urv is formally obtained from V^surf by replacing the surface 
tension a by an effective term ^/{2R). 



Appendix B 

In this appendix we briefly present the results of Ref. |33] for the level density in the sphere and we 
make explicit the different terms appearing in Eq. (^ for the perfect billiard. 

The periodic orbits in the sphere are regular polygons in diametral planes. They are labeled by 
two numbers (n,t), n being the number of sides and t the winding number of the orbit around the 
center (n > 2t). Note that n is here the same as the number of bouncing points appearing in the main 
text (Eq. (^)). The oscillating part of the level density in the sphere reads : 

+ 00 +00 

Posc{k) = 51 51 ■^n,t{k) sm{kLn,t + t'n,t7r/2) . (Bl) 

4=1 n=2t 

The shortest orbits are the pendulating orbit (n = 2, t = 1), the triangle (n = 3, t = 1), and 
the square (n = 4, t = 1). The triangle and the square are sufficient to understand the qualitative 



features of the shell and supershell structure (see e.g. |5C,5l[). Each orbit bounces on the surface with 
a constant normal angle 9n,t = (1 ~ 2t/n)7r/2 and has a length L„^f = 2nRcos9n,t- The pendulating 
orbit occurs in a two-parameters family (the parameters determine the direction of bouncing) whereas 
all the other orbits form three-parameters families. Hence the bouncing ball (with n = 2t) has to be 
treated separately. The explicit formulae for A{k) and u in Eqs. (^) and ( |BlD are {t>l) : 




if n = 2t , 
if n>2t , 



(B2) 



-^^^ if n = 2t, 



and 



2ds{-lY s\n{2en,t)\j^^^R{Rkf/^ if n > 2t , 

where = 2 is the spin degeneracy. 

One sees that the amplitude corresponding to the pendulating orbit is proportional to k whereas 
the other families have a larger weight (proportional to k'^/'^). Generally speaking, one can show that 
the contribution of a d-parameter family has an extra power with respect to that of an isolated 
orbit (51 . 
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Appendix C 



In this appendix we derive approximate analytical expressions for the oscillatory part of the total 
energy and of the free energy at finite temperature. Similar results concerning the entropy, the free 
energy ete. . . have been previously obtained by Kolomiets, Magner and Strutinsky in Refs. 
We nevertheless briefly outline the derivation of the formulae because the references just quoted are 
not very explicit and difficult to follow. The formulae are derived in the framework of a general PO 
expansion : the level density is noted p{e), it is separated in a smooth term p{e) and an oscillating term 
Posc(e)- In the present work, we denote all the smooth terms with an upper bar and the oscillating 
terms with a subscript "osc", except for the oscillating part of the energies which have a subscript 
"shell" according to the general convention in the field. /9osc(e) is supposed to be of the form : 

PoUe)=ReJ2B{e)e'^('y^ , (CI) 

PO 

where S{e) is the action of the PO considered. B{e) is an orbit dependent amplitude which is of 
order for chaotic systems, of order h~'^ for typical integrable systems in three dimensions and of 
order U^^^"^ for rotationally symmetric systems as the spherical billiard where families of orbits are 
characterized by 3 parameters [^2| . 

We will estimate the asymptotic form of several integrals, all of the same type, and we first display 
a formula often used below. Let g{e) be a slowly varying function of e (as B{e) is supposed to be), 
and g' its first derivative. Let 0(e — fi) = [1 + exp{(e — fi)/kgT}]~^ be the Fermi function, fj. being the 
chemical potential. If S{p) ^ h one has : 



p., - ...... _ I ^ _ ^ i;;^,,., , 

(C2) 

where the integral has been evaluated by a contour integration in the complex plane (see e.g. |[55[|). 
In the evaluation of the integral we have neglected the contribution of a part of the contour located 
on the positive imaginary axis, this is legitimate provided the temperature is small compared to the 
Fermi energy (degenerate Fermi gas approximation). X = TTS'{p)kgT/h is a dimensionless quantity 
which can be considered as small if the period 5' of the orbit is small compared to a characteristic 
thermal time h/kgT. Fi, F2 and F-^ are dimensionless damping functions : 



, , X , , X^coshX , , , sinh^X, , , 



For obtaining the total energy starting from (CI), one first determines the chemical potential /i 



through the equality N = AA(/i), where is the number of electrons and M{fi) = J p{e)4){e — p)de. 
As /9(e), M can be separated in a smooth term M (the Weyl term) plus an oscillating part AAqsc- 
Accordingly, p, can be separated in a smooth function of plus an oscillating term : p = p, + Posc, 
where A'' = M{p). 

Then the total electronic energy is E = J e0(e — p)p{e)de. It can also be separated into a smooth 
part E and an oscillating part which is denoted -Esheii throughout the paper (more precisely £'gheii(-^) T) 
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in presence of temperature) in accordance with the general notations in the field. -Egheii reads appro- 
ximatively : 



EsheU = J e p{e)(j){e - n)de - J e p{e)4i{e - fL)d€ 

- J e ['/'(e - fi) - /^osc 0'(e - fl)]de - J e p(e)0(e - fl)de (C4) 
= J e posc{e)(t>{e - fJ.)de + Pose (^^^ _- Pose J {e - p)p{e)(l)' {e - p)de . 

The last term of the r.h.s. of ( p4D is subdominant, moreover it is zero at zero temperature, hence 



we drop it in the following. Then, from (CI) and (|C2| ) one obtains : 



i?sheii(iV,T) ^ -Re ^ i^-A-^'B{fl)F2{X)e'^(f^y^ , (C5) 
where X is computed as X with p replacing p. 

The free energy F(N, T) is a quantity more appropriate than the total energy to evaluate the 
effects of electronic temperature on the abundance of clusters in the beam (see the discussion in the 
main text). It is defined by 

F{N,T) = pN + £°°dep{e)^{e-p) where ^{e - p) = -kBT\xi[l + e^^^ - ^) I^^T^ . (C6) 
The oscillating part of the free energy is denoted by Fsheii(-/V, T) and can be evaluated similarly to 



what has been done in (C4). This yields : 



Fshcii(iV,T) ^ - Re ^ ( B{-p)F,{X)e^Sm/^ , (C7) 



PO 



In the particular case of a billiard whose level density is of the type (|B1|), Eq. (P?]) reads : 



h'^kl „ 2Aik„) 



where k^ is defined hy p = fi'^k'j^/ {2m). A formula of this type seems to have been derived first by 
Dingle in Ref. |56|. We have verified that this formula is of very good accuracy in the spherical billiard 



(see Fig. 4). An equally good agreement is obtained for the comparison of -Esheii (as given by (p5|)) 
with the exact result. For relatively low values of N^/^ (say N^/'^ < 6), an even better agreement 
can be obtained by still using the semi-classical level density, but evaluating integrals such as (|C^ ) 
numerically. 

In the sphere, the main contribution to ( |C8| ) comes from orbits occurring in three-parameters 
families (c/ Appendix B). Considering that k^ is of order Kp ~ l/r^ and that L and R scale like 
rs N^/^ one obtains the following leading order : A{kn) ~ R{Rkp)^/'^ ~ VgN^/^. Hence the schematic 



large behaviour of (|C8D reads 
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Figure 4: Fsheii (expressed in units of the bulk Fermi energy Sp) as a function of N'^^^ in a perfectly 
spherical aggregate. FghcU is denoted Fshcii{^,T) in the main text. The different plots correspond 
to temperatures r = 0.01 (upper plot), r = 0.02 and 0.03 (lower plots), where r = kgT/ep. For 
sodium these values correspond to T = 376 K, 752 K and 1128 K. The solid lines correspond to the 
determination of -Fgheii obtained by using the exact spectrum of the spherical billiard, and the dashed 
lines to Eq. 



shell 



epN'^/^Fi{TN'^/^)sm{N^/^) . (C9) 



where r = kgT/ep is the temperature expressed in units of the bulk Fermi energy. Here we have 
dropped for clarity important but dimensionless factors in the sine and -Fi function : we just want to 
illustrate the typical N dependence of KjhcU- We have adopted the same type of notation in the text 
(^^- The behaviour ( |C9| ) is in agreement with the findings of Kolomiets, Magner and Strutinsky in 
Refs. [^,54 1 and with Ref. |jl8| where Bohr and Mottelson used schematic forms of the level density. 



We emphasize that we have used here a generic PO expansion, and Eqs. (C5,^) are valid for any 



system (chaotic or integrable) of independent fermions moving in an external potential. 

Note finally that we have here computed the thermodynamical quantities in the grand canonical 
ensemble. The number of electrons in a cluster being exactly conserved, the canonical description 
should be used instead (hence we should have noted -F(/x, T) instead of F{N, T) : in all the appendix 
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should be understood as the mean number of electrons). The difference between the two ensembles 
has been studied in Ref. where it is shown to give discrepancies of order of 0.05 eV (or 0.1 eV at 
best) in the free energy difference F{N — 1) — F{N). This difference is expected to decrease in the 
large limit and moreover it plays no role in the discussion of the effects of temperature given in the 
main text. 
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